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Abstract 

In  this  paper,  a  three-dimensional  PEM  fuel  cell  model  with  a  consistent  water  transport  treatment  in  the  membrane  electrode  assembly  (MEA) 
has  been  developed.  In  this  new  PEM  fuel  cell  model,  the  conservation  equation  of  the  water  concentration  is  solved  in  the  gas  channels,  gas  diffusion 
layers,  and  catalyst  layers  while  a  conservation  equation  of  the  water  content  is  established  in  the  membrane.  These  two  equations  are  connected 
using  a  set  of  internal  boundary  conditions  based  on  the  thermodynamic  phase  equilibrium  and  flux  equality  at  the  interface  of  the  membrane 
and  the  catalyst  layer.  The  existing  fictitious  water  concentration  treatment,  which  assumes  thermodynamic  phase  equilibrium  between  the  water 
content  in  the  membrane  phase  and  the  water  concentration,  is  applied  in  the  two  catalyst  layers  to  consider  water  transport  in  the  membrane  phase. 
Since  all  the  other  conservation  equations  are  still  developed  and  solved  in  the  single-domain  framework  without  resort  to  interfacial  boundary 
conditions,  the  present  new  PEM  fuel  cell  model  is  termed  as  a  mixed-domain  method.  Results  from  this  mixed-domain  approach  have  been 
compared  extensively  with  those  from  the  single-domain  method,  showing  good  accuracy  in  terms  of  not  only  cell  performances  and  current 
distributions  but  also  water  content  variations  in  the  membrane. 

©  2006  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Many  multi-dimensional  PEM  fuel  cell  models  have  been 
developed  in  the  past  decade  to  facilitate  cell  design  and  opti¬ 
mization.  These  models  generally  fall  into  two  categories:  multi- 
domain  and  single-domain  method.  The  studies  of  Gurau  et 
al.  [1]  and  Berning  et  al.  [2]  were  based  on  the  multi-domain 
method,  in  which  the  computational  domain  was  divided  into 
a  number  of  sub-domains  and  different  sets  of  conservation 
equations  were  developed  in  different  sub-domains.  Interfacial 
boundary  conditions  were  further  established  to  connect  these 
equations.  In  the  single-domain  method,  one  set  of  conservation 
equations  was  applied  to  different  regions  of  a  PEM  fuel  cell. 
In  order  to  switch  on/off  a  specific  equation  in  a  specific  region, 
special  numerical  treatments  have  been  used,  including  defining 
extremely  large  or  small  physical  and  transport  parameters  in 
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the  region,  i.e.  very  small  proton  conductivity  outside  the  mem¬ 
brane  electrode  assembly  (MEA).  The  studies  of  Um  et  al.  [3], 
Dutta  et  al.  [4],  Siegel  et  al.  [5],  and  Mazumder  and  Cole  [6] 
were  in  this  category. 

Since  there  is  only  one  set  of  conservation  equations  and 
no  interfacial  boundary  condition,  the  single-domain  method  is 
easier  to  formulate  and  implement  into  a  general  CFD  package. 
Recently,  this  method  has  experienced  rapid  development.  For 
example,  it  has  been  applied  to  study  electron  transport  phe¬ 
nomena  [7,8],  heat  transfer  [9],  large-scale  simulations  [10-12], 
and  two-phase  flows  and  dynamics  [13,14].  Extensive  model 
validations  have  also  been  conducted  [15,16]. 

A  weakness  of  the  single-domain  method  lies  in  its  inability  to 
handle  water  transport  through  the  membrane  phase  directly,  as 
theoretically  the  water  content  has  to  be  solved  in  the  membrane 
phase  in  MEA  while  the  water  concentration  has  to  be  solved  in 
the  other  regions.  In  the  study  of  Dutta  et  al.  [4],  the  MEA  region 
was  completely  neglected  from  the  computational  domain.  As 
such,  a  simplified  treatment  was  applied  for  water  transport  in 
the  membrane  phase.  In  the  work  of  Um  et  al.  [3],  a  fictitious 
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Nomenclature 

a  water  activity 

c  molar  concentration  (mol  m-3) 

Cp  constant-pressure  heat  capacity  (J  (kg  K)  - 1 ) 

D  mass  diffusivity  (m2  s-1) 

D\  water  content  diffusivity  (mol (ms)-1) 

EW  equivalent  weight  of  the  membrane  (kg  mol-1 ) 

F  Faraday  constant  (96487  C  mol- 1 ) 

i  current  density  vector  (AM-2) 

4vg  average  current  density  (A  m-2) 

k  thermal  conductivity  (W  (m  K)- 1 ) 

rid  electro-osmotic  drag  coefficient 

Aw  net  water  flux  through  the  membrane 

(mol  (m2  s)-1) 
p  pressure  (Pa) 

q  interfacial  flux 

Ru  universal  gas  constant  (J  (mol  K)  - 1 ) 

S  source  term 

T  temperature  (K) 

u  velocity  (ms-1) 

Greek  symbols 

a  net  water  transfer  coefficient 

s  porosity 

sm  fraction  of  the  membrane  phase  in  the  catalyst 
layer 

0  phase  potential  (V) 

k  proton  conductivity  (S  m- 1 ) 

A  water  content 

/x  chemical  potential 

p  density  (kg  m- 3 ) 

o  electronic  conductivity  (S  m-1) 

r  viscous  stress  tensor 

Supercripts 
cl  catalyst  layer 

eff  effective  value 

m  membrane 

Subscripts 

cl  catalyst  layer 

e  electrolyte  or  energy 

g  gaseous  phase 

i  species 

m  membrane 

s  electron 

sat  saturation  value 

w  water 


water  concentration  was  derived  to  replace  the  variable  of  water 
content  and  therefore  a  same  form  of  water  transport  equation 
could  be  solved  in  a  single-domain  framework  [17].  The  same 
treatment  has  also  been  proposed  by  Kulikovsky  [18],  but  it  was 
applied  only  in  the  two  catalyst  layers.  Since  in  this  treatment, 
the  water  content  in  the  membrane  phase  was  assumed  to  be  in 


the  thermodynamic  equilibrium  with  the  water  concentration  in 
every  location,  including  inside  the  membrane,  it  seems  that  it 
is  more  appropriate  to  apply  this  treatment  only  in  the  catalyst 
layers,  as  the  membrane  phase  is  pervasively  distributed  and  in 
extensive  contact  with  water  in  these  two  regions  [18].  It  lacks 
theoretical  basis  to  apply  this  method  inside  the  membrane  and 
this  treatment  is  thus  an  approximation  in  this  region.  However, 
this  treatment  has  been  widely  used  in  the  studies  of  Meng  and 
Wang  [7,8,10],  Ju  et  al.  [9],  and  Wang  and  Wang  [12,19],  and 
it  has  been  successfully  validated  by  Ju  and  Wang  [15]  and  Ju 
et  al.  [16]  in  terms  of  cell  performances  and  current  distribu¬ 
tions.  In  the  study  of  Siegel  et  al.  [5],  an  extra  transport  equation 
for  the  dissolved  water  concentration  was  established  in  MEA, 
and  a  source  term  in  the  form  of  convective  mass  transfer  was 
used  to  account  for  the  water  dissolution  rate  into  the  membrane 
phase.  However,  no  expression  has  been  presented  for  calculat¬ 
ing  this  parameter.  It  is  not  clear  how  water  transport  through  the 
membrane  was  handled  in  the  study  of  Mazumder  and  Cole  [6] . 

In  this  paper,  a  three-dimensional  PEM  fuel  cell  model  with 
a  consistent  water  transport  treatment  in  the  membrane  elec¬ 
trode  assembly  has  been  developed.  The  concept  of  the  fictitious 
water  concentration  is  applied  in  the  two  catalyst  layers.  An 
extra  conservation  equation  of  the  water  content  is  solved  in  the 
membrane.  Unlike  Kulikovsky  [18],  however,  a  different  set  of 
internal  boundary  conditions  at  the  interface  of  the  membrane 
and  the  catalyst  layer  has  been  established,  extending  the  one¬ 
dimensional  formulation  of  Springer  et  al.  [20]  into  the  present 
three-dimensional  framework.  Since  all  the  other  conservation 
equations  are  still  developed  and  solved  in  the  single-domain 
paradigm,  the  present  PEM  fuel  cell  model  is  termed  as  a  mixed- 
domain  approach.  In  this  paper,  results  from  this  mixed-domain 
approach  have  been  compared  to  those  from  the  single-domain 
method  in  detail. 

2.  Theoretical  formulation 

The  conservation  equations  of  mass,  momentum,  species 
concentration,  proton,  electron,  and  energy  are  still  formulated 
exactly  in  the  same  forms  as  in  the  single-domain  framework. 
They  are  in  the  following  forms: 


Mass: 

V  •  (pii)  =  0 

(1) 

Momentum: 

1 

-^-V  •  ( puu )  =  —S/p  +  V  •  r  +  Su 

£z 

(2) 

Species: 

V  •  (Sc,-)  =  V  •  (DfVa)  +  Si 

(3) 

Proton  transport: 

V  •  (/CeffV0e)  +  Se  =  0 

(4) 

Electron  transport: 

V  •  (CTeffV0s)  +  Ss  =  0 


(5) 
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Energy: 

V-(pCpuT)  =  V  •  (feeffvr)  +  5r  (6) 


In  these  equations,  the  source  terms  and  the  other  relevant 
physicochemical  parameters  are  presented  in  Meng  and  Wang 
[7,10]  and  Ju  et  al.  [9]. 

Among  the  species  concentration  equations,  the  equation  of 
the  hydrogen  concentration  is  solved  only  on  the  anode  side 
while  the  equation  of  the  oxygen  concentration  on  the  cathode 
side.  In  the  present  mixed-domain  method,  the  equation  of  the 
water  concentration  is  solved  only  in  the  gas  channels,  gas  diffu¬ 
sion  layers  (GDL),  and  catalyst  layers  on  both  anode  and  cathode 
sides.  Furthermore,  in  the  two  catalyst  layers,  the  fictitious  water 
concentration  has  been  included  to  take  into  account  of  water 
transport  in  the  membrane  phase  [17,18].  The  effective  water 
diffusion  coefficient  in  the  two  catalyst  layers  can  be  expressed 
as 


ncl.eff  =  1.5  Dcl,g  ,  1.5  D  RuT^X 


P sat  da 


More  details  can  be  found  in  Um  and  Wang  [17]  and  Kulikovsky 
[18]. 

Unlike  the  prior  single-domain  method,  a  conservation  equa¬ 
tion  of  the  water  content  is  solved  in  the  membrane  in  this 
mixed-domain  method.  The  conservation  equation  is  derived 
as  follows: 


V-(D1VA)  +  ^  =  0 


where  the  source  term  arising  from  the  electro-osmotic  drag  can 
be  expressed  as 


In  Eq.  (8),  because  the  fluid  velocity  is  very  small  inside  the 
membrane,  the  convective  effect  has  been  neglected. 

As  in  the  prior  studies  of  Meng  and  Wang  [7,8,10]  and  Um 
and  Wang  [17],  the  water  content  diffusivity  in  Eqs.  (7)  and  (8) 
can  be  estimated  as 


Pm  _  Pm 
EW  w  EW 

'  3.1  x  10-7A(ea28A  -  1)  •  e[-2346/rl  0  <  A  <  3 

4.17  x  10_8A(1  +  161e-A)  •  e^-2346/7^  otherwise 

(10) 


Its  variation  is  shown  in  Fig.  1. 

The  conservation  equation  of  the  water  content,  Eq.  (8),  is 
connected  to  the  water  concentration  equation  by  a  set  of  internal 
boundary  conditions  at  the  two  interfaces  between  the  mem¬ 
brane  and  the  catalyst  layers  on  both  anode  and  cathode  sides. 
As  shown  in  Fig.  2,  the  interfacial  boundary  conditions  can  be 
established  based  on  the  thermodynamic  phase  equilibrium  and 
the  flux  equality.  The  general  thermodynamic  phase  equilibrium 
conditions  at  the  interfaces  are  [21] 

Tm  =  Tcl  (11a) 


Fig.  1.  Variation  of  water  content  diffusivity  in  the  membrane  phase. 


Pm  =  Pd 


Mw  =  Mw 


(Hb) 

(11c) 


The  general  flux  equality  conditions  include  the  equal  energy 
and  water  fluxes  at  the  interface 


qf  =  q?  (12a) 

<7w  =  7w  (12b) 

Since  the  single-domain  method  is  used  for  solving  the  energy 
equation,  Eqs.  (11a)  and  (12a)  will  be  closely  approached  but 
are  not  explicitly  required.  Eq.  (lib)  is  required  for  solving 
water  transport  in  the  membrane  once  liquid  water  is  involved, 
as  suggested  in  Weber  and  Newman  [22].  In  the  present  model 
development  and  the  following  numerical  calculations,  we  will 


MEM 


i 


Fig.  2.  Schematic  of  an  interface  and  the  interfacial  boundary  conditions. 
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only  consider  the  pseudo  single-phase  cases,  and  therefore,  this 
equation  can  be  neglected. 

According  to  Springer  et  al.  [20],  Eq.  (11c)  can  be  approxi¬ 
mately  replaced  using  the  following  empirical  expression: 

f  0.043  +  17.18a  —  39.85 a2  +  36.0a3  0  <  a  <  1 

A  =  <  “  (13) 

[14+1.4(a-l)  1  <  a  <  3 

where  the  parameter  X  represents  the  water  content  on  the  mem¬ 
brane  side  of  the  interface  and  the  parameter  a  is  the  water 
activity  on  the  catalyst  layer  side  of  the  interface. 


The  water  flux  equality  at  the  interface  can  be  further  derived 
as 

\-DkVX  + 

Eqs.  (13)  and  (14)  constitute  a  set  of  the  interfacial  boundary 
conditions  for  connecting  the  two  conservation  equations  of  the 
water  content  in  the  membrane  and  the  water  concentration  in 
the  other  regions. 

The  above  conservation  equations,  Eqs.  (l)-(6)  and  (8),  and 
the  set  of  interfacial  boundary  conditions,  Eqs.  (13)  and  (14), 
complete  the  model  development  in  the  present  mixed-domain 


ndi 

F 


=  -D 


cl,eff 


w 


vcw  + 


ndi 

F 


m 


(14) 


cl 


1 3400.0 

13027.3 

12654.5 
12281.8 

11909.1 

11536.4 

11163.6 
10790.9 

10418.2 

10045.5 
9672.7 
9300.0 


y 


13400.0 

13027.3 

12654.5 
12281.8 

11909.1 

11536.4 

11163.6 
10790.9 

10418.2 

10045.5 
9672.7 
9300.0 


Fig.  3.  Current  distribution  in  the  mid-thickness  of  the  membrane:  (a)  from  the  mixed  domain  method,  (b)  from  the  single-domain  method  (unit:  Am  2). 
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method.  The  other  relevant  physicochemical  relationships, 
boundary  conditions,  and  numerical  treatments  follow  the  same 
procedures  as  in  the  single-domain  approach,  as  detailed  in 
Refs.  [3,7,10].  This  numerical  model  has  been  implemented 
into  a  general  CFD  package,  Fluent.  The  computational  grid 
used  in  the  present  study  is  the  same  as  in  the  previous  work 
[7,8],  which  is  generated  based  on  careful  grid  independence 
studies.  In  the  next  section,  results  from  this  mixed-domain 
approach  will  be  compared  with  those  from  the  single-domain 
method,  which  have  been  extensively  validated  in  terms  of  cell 
performances  and  current  distributions  [15,16]. 

3.  Result  and  discussion 

The  present  numerical  calculations  are  conducted  using  a 
single  straight-channel  PEM  fuel  cell  with  a  thin  membrane  of 
25  pm.  Its  geometry  and  the  other  related  parameters  are  pre- 


tion. 


Table  1 

Inlet  humidification  temperature  and  relative  humidity  at  cell  temperature  of 
80  °C 


Case 

number 

Anode 

Cathode 

Humidification 
temperature  (°C) 

Relative 

humidity 

(%) 

Humidification 
temperature  (°C) 

Relative 

humidity 

(%) 

Case  1 

80 

100 

20 

5 

Case  2 

50 

26 

50 

26 

Case  3 

80 

100 

80 

100 

sented  in  Meng  and  Wang  [7,8].  In  the  present  configuration,  the 
T-coordinate  is  in  the  through-membrane  direction,  y-coordinate 
the  along-channel  direction,  z-coordinate  the  lateral  direction. 
In  order  to  make  a  comprehensive  comparison,  three  different 
cases  are  designed  for  the  present  numerical  study,  as  shown 
in  Table  1 .  Since  the  cell  operates  at  a  constant  temperature  of 
80  °C,  two  of  the  cases  are  in  low-humidity  operation  conditions, 
consistent  with  the  present  pseudo  single-phase  calculations. 


Fig.  6.  Variation  of  water  content  inside  the  membrane  under  the  gas  channel: 
(a)  at  the  inlet  region,  (b)  at  the  outlet  region. 
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The  fully  humidified  case  is  included  for  completeness.  The 
single-domain  calculations  are  based  on  the  model  presented  in 
Meng  and  Wang  [7,10],  which  applies  the  fictitious  water  con¬ 
centration  concept  in  ME  A  [17].  In  both  the  mixed-domain  and 
the  single-domain  methods,  the  anisotropic  electron  transport 
phenomenon  is  handled  using  the  simplified  method  developed 
in  Meng  [23],  but  the  contact  resistance  is  neglected. 

Fig.  3  illustrates  current  distributions  in  the  mid-thickness 
of  the  membrane  from  both  the  present  mixed-domain  and  the 
single-domain  methods  for  case  1.  Results  show  an  excellent 
agreement  at  the  first  two-third  length  of  the  cell  from  the  cell 
inlet.  The  current  distribution  from  the  single-domain  model 
shows  irregular  variation  starting  from  the  two-third  cell  length. 
The  calculated  current  distribution  from  the  present  mixed- 
domain  method  gives  smooth  variation  in  the  entire  region. 
Variations  of  the  average  current  density  in  the  along-channel 
direction  from  the  two  methods  are  further  compared  in  Fig.  4. 
The  curves  in  Fig.  4  vary  in  the  exactly  same  trend  and  show 
a  very  good  agreement.  In  consistency  with  Fig.  3,  the  aver¬ 
age  current  density  from  the  single-domain  method  shows  slight 
oscillations  starting  from  the  two-third  cell  length.  The  reason 
for  the  oscillations  will  be  explained  later  in  this  section. 

Fig.  5  presents  variations  of  the  amount  of  water  transported 
through  the  membrane  in  terms  of  a  net  water  transfer  coefficient 
a,  which  is  defined  as  the  ratio  of  the  local  net  water  transfer  rate 


Fig.  7.  Variation  of  the  derivative  term,  dX/da,  with  water  content. 


through  the  membrane  and  the  average  water  production  rate  in 
the  cathode  catalyst  layer 


a  = 


2FNW 

7avg 


(15) 


Fig.  8.  Water  content  distribution  in  a  cross  section  perpendicular  to  the  membrane:  (a)  at  the  inlet  region,  (b)  in  the  middle,  (c)  at  the  outlet  region. 
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In  Eq.  (15),  a  positive  value  means  that  water  is  transferred  from 
the  anode  to  the  cathode  side.  The  net  water  transfer  through  the 
membrane  is  caused  by  two  mechanisms  in  the  present  calcu¬ 
lations,  namely  the  electro-osmotic  drag  from  the  anode  to  the 
cathode  side  and  water  diffusion,  whose  direction  depends  on 
the  gradient  of  the  water  concentrations  on  the  two  sides. 

In  Fig.  5,  variations  of  the  net  water  transfer  coefficient  calcu¬ 
lated  from  the  two  methods  are  in  an  excellent  agreement  before 
the  result  from  the  single-domain  method  starts  to  oscillate.  In 
both  methods,  the  calculated  net  water  transfer  is  from  the  anode 
to  the  cathode  side  at  the  beginning  of  the  cell  because  of  the 
dual  effects  of  the  higher  water  concentration  on  the  anode  side, 


Fig.  9.  Water  content  distribution  in  the  mid-thickness  of  the  membrane. 


which  results  in  the  forward  water  diffusion,  and  the  electro- 
osmotic  drag.  At  around  y  =  4  cm,  the  net  water  transfer  coeffi¬ 
cient  changes  to  negative  values.  This  is  dictated  by  the  backward 
water  diffusion  caused  by  the  higher  water  concentration  on  the 
cathode  side,  which  is  produced  by  the  electrochemical  reaction. 

Fig.  6  shows  water  content  variations  inside  the  membrane 
under  the  middle  of  the  gas  channel.  In  Fig.  6a  close  to  the  inlet 
region,  the  water  content  has  higher  values  on  the  anode  side, 
resulting  in  the  forward  water  diffusion.  Results  from  the  two 
methods  are  in  an  excellent  agreement.  In  Fig.  6b  close  to  the 
outlet  region  of  the  cell,  the  water  content  shows  higher  val¬ 
ues  on  the  cathode  side,  resulting  in  the  strong  backward  water 
diffusion.  An  interesting  point  is  that  although  the  curves  show 
nearly  the  same  gradient  in  the  majority  of  the  membrane,  the 
result  from  the  single-domain  approach  shows  a  sharp  decrease 
at  a  water  content  value  of  X  =  14  close  to  the  cathode  side.  This 
is  not  a  physical  phenomenon  but  a  numerical  result  caused  by 
the  fictitious  water  concentration  treatment,  in  which  the  water 
diffusivity  involves  a  derivative  term  in  the  form  of  df/da.  Fig.  7 


Fig.  10.  Variation  of  water  concentration  in  the  gas  channel:  (a)  on  the  anode 
side,  (b)  on  the  cathode  side. 
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provides  the  variation  of  this  derivative  term  calculated  using  Eq. 
(10).  It  shows  an  abrupt  decrease  at  a  water  content  of  14.  This 
result  indicates  that  the  fictitious  water  concentration  approxi¬ 
mation  could  cause  an  incorrect  water  content  variation  inside 
the  membrane.  In  fact,  the  sharp  decrease  of  the  water  con¬ 
tent  in  the  single-domain  approach,  as  shown  in  Fig.  6b,  results 
in  the  lower  water  content  values  in  the  majority  of  the  mem¬ 
brane,  which  in  turn  causes  the  lower  water  content  diffusivity, 
as  can  be  seen  in  Fig.  1,  and  consequently  the  higher  net  water 
transfer  coefficient  at  the  outlet  region  in  Fig.  5.  In  addition, 
this  incorrect  water  content  variation  also  causes  oscillations 
in  the  current  distribution  and  the  net  water  transfer  coefficient 
in  Figs.  3-5.  It  should  be  noted  that  in  this  case  this  type  of 
problem  only  occurs  under  the  fully  humidified  condition  with 
a  water  activity  above  unity,  i.e.  close  to  the  outlet  region,  but 
the  incorrect  water  content  variation  could  also  occurs  under 
low-humidification  conditions,  as  discussed  later  in  this  section. 

Water  content  distributions  calculated  from  the  present 
mixed-domain  model  are  illustrated  in  Figs.  8  and  9.  In  Fig.  8, 


(b)  Channel  Length  (cm) 


Fig.  11.  Variation  of  the  average  current  density  in  the  along-channel  direction: 
(a)  from  case  2,  (b)  from  case  3. 


water  content  distributions  are  shown  in  three  different  cross 
sections  perpendicular  to  the  membrane.  These  results  clearly 
show  that  the  water  content  distribution  in  the  membrane  varies 
not  only  from  the  anode  to  the  cathode  side  but  also  in  the  lateral 
direction,  with  the  higher  water  content  under  the  land  than  that 
under  the  channel.  Fig.  8a  illustrates  that  at  the  inlet  region  of 
the  cell,  although  the  water  content  is  higher  on  the  anode  side 
under  the  gas  channel,  as  also  presented  in  Fig.  6a,  it  is  higher  on 
the  cathode  side  under  the  land.  Fig.  9  presents  the  water  content 
distribution  in  the  mid-thickness  of  the  membrane,  which  clearly 
shows  that  the  water  content  increases  from  the  inlet  to  the  outlet 
region,  painting  a  complete  picture  of  the  water  content  variation 
inside  the  membrane. 

Average  water  concentration  variations  in  both  the  anode 
and  the  cathode  gas  channels  from  the  two  methods  are  pre¬ 
sented  in  Fig.  10.  Results  from  the  two  methods  are  in  very  good 
agreements  as  expected.  Results  are  also  consistent  with  the  net 
water  transfer  coefficients  presented  in  Fig.  5.  For  example,  the 
water  concentration  in  the  anode  gas  channel  initially  decreases 
owing  to  the  dual  effects  of  the  forward  water  diffusion  and  the 


Fig.  12.  Variation  of  the  net  water  transfer  coefficient:  (a)  from  case  2,  (b)  from 
case  3. 
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electro-osmotic  drag,  and  then  increases  dominated  by  the  back¬ 
ward  water  diffusion.  Furthermore,  the  differences  from  the  two 
methods  are  consistent  on  the  anode  and  the  cathode  sides  as 
well. 

Variations  of  the  average  current  density  from  the  mixed- 
domain  and  single-domain  methods  are  further  compared  in 
Fig.  11,  showing  very  good  agreements  for  both  cases  2  and  3. 
The  net  water  transfer  coefficient  from  the  two  methods  is  com¬ 
pared  in  Fig.  12.  Results  from  case  2  in  Fig.  12a  show  a  very 
good  agreement.  Although  results  from  case  3  in  Fig.  12b  are  in 
the  same  trend,  the  result  from  the  single-domain  method  shows 
large  oscillations.  These  comparisons  indicate  that  the  single¬ 
domain  method  performs  better  under  the  low-humidification 
conditions  than  under  the  fully  humidified  ones,  in  consis¬ 
tency  with  the  result  from  case  1 .  However,  even  under  a  low- 
humidification  condition  in  case  2,  the  calculated  water  content 
in  the  membrane  from  the  single-domain  approach  is  different 
from  the  mixed-domain  method,  as  shown  in  Fig.  13a  at  a  loca¬ 
tion  directly  under  the  gas  channel  close  to  the  outlet  end.  This 
can  also  be  attributed  to  the  approximate  fictitious  water  concen- 


Fig.  13.  Variation  of  water  content  inside  the  membrane  under  the  gas  channel 
at  the  outlet  region:  (a)  from  case  2,  (b)  from  case  3. 


tration  treatment  in  the  membrane.  Under  the  fully  humidified 
condition  in  case  3,  a  rapid  decrease  of  the  water  content  cal¬ 
culated  from  the  single-domain  method  again  occurs  inside  the 
membrane,  as  shown  in  Fig.  13b. 

Based  on  a  comprehensive  numerical  analysis  presented 
herein,  it  seems  that  the  single-domain  method  using  the  fic¬ 
titious  water  concentration  approximation  should  be  improved 
as  it  lacks  theoretical  basis  inside  the  membrane  and  gives  incor¬ 
rect  water  content  variations,  although  it  can  provide  correct  cell 
performances  and  current  variations.  The  present  mixed-domain 
method,  which  applies  a  set  of  interfacial  boundary  conditions 
to  connect  the  conservation  equation  of  the  water  content  in  the 
membrane  and  the  water  concentration  equation  in  the  other 
regions,  is  established  on  a  solid  theoretical  basis  and  provides 
correct  results  in  terms  of  not  only  cell  performances  and  current 
variations  but  also  water  content  variations  inside  the  membrane. 
Furthermore,  it  seems  that  model  validations  based  on  cell  per¬ 
formances  and  current  variations  are  still  insufficient  since  small 
differences  in  the  current  variation  between  the  numerical  and 
experimental  data  could  indicate  an  incorrect  result  but  it  could 
easily  escape  a  researcher’s  attention.  It  seems  that  the  water 
content  in  the  membrane  and  its  related  parameters,  i.e.  the 
membrane  resistance,  might  have  to  be  used  for  better  model 
validation. 

4.  Conclusion 

In  this  paper,  a  three-dimensional  PEM  fuel  cell  model  with 
a  consistent  water  transport  treatment  in  the  membrane  elec¬ 
trode  assembly  has  been  developed.  In  this  new  PEM  fuel  cell 
model,  the  conservation  equation  of  the  water  concentration  is 
solved  only  in  the  gas  channels,  gas  diffusion  layers,  and  cata¬ 
lyst  layers  while  a  conservation  equation  of  the  water  content  is 
established  in  the  membrane.  These  two  equations  are  connected 
using  a  set  of  interfacial  boundary  conditions  based  on  the  ther¬ 
modynamic  phase  equilibrium  and  flux  equality  at  the  interface 
of  the  membrane  and  the  catalyst  layer.  In  fact,  this  theoretical 
formulation  extends  the  Springer’s  one-dimensional  interfacial 
treatment  into  a  three-dimensional  framework.  In  addition,  the 
prior  fictitious  water  concentration  treatment  is  applied  in  the 
two  catalyst  layers  on  a  solid  theoretical  basis  to  consider  water 
transport  in  the  membrane  phase.  Since  all  the  other  conser¬ 
vation  equations  of  mass,  momentum,  species  concentration, 
proton  and  electron  transport,  and  energy  are  still  solved  in  the 
single-domain  framework  without  resort  to  interfacial  bound¬ 
ary  conditions,  the  present  new  PEM  fuel  cell  model  is  termed 
as  a  mixed-domain  method.  Results  from  this  mixed-domain 
method  have  been  compared  extensively  with  those  from  the 
single-domain  approach.  The  present  model  provides  accurate 
numerical  results  in  terms  of  not  only  cell  performances  and 
current  variations  but  also  water  content  variations  inside  the 
membrane. 
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